CONSISTENT NEWTON-RAPHSON VS. FIXED-POINT FOR VARIATIONAL 
MULTISCALE FORMULATIONS FOR INCOMPRESSIBLE NAVIER STOKES 



D. Z. TURNER, K. B. NAKSHATRALA, AND K. D. HJELMSTAD 

Abstract. The following paper compares a consistent Newton- Raphson and fixed-point itera- 
tion based solution strategy for a variational multiscale finite element formulation for incompress- 
ible Navier-Stokes. The main contributions of this work include a consistent linearization of the 
Navier-Stokes equations, which provides an avenue for advanced algorithms that require origins 
in a consistent method. We also present a comparison between formulations that differ only in 
their linearization, but maintain all other equivalences. Using the variational multiscale concept, 
we construct a stabilized formulation (that may be considered an extension of the MINI element to 
nonlinear Navier-Stokes). We then linearize the problem using fixed-point iteration and by deriving 
a consistent tangent matrix for the update equation to obtain the solution via Newton-Raphson 
iterations. We show that the consistent formulation converges in fewer iterations, as expected, 
for several test problems. We also show that the consistent formulation converges for problems 
for which fixed-point iteration diverges. We present the results of both methods for problems of 
Reynold's number up to 5000. 



1. INTRODUCTION 

Within the context of variational multiscale formulations for incompressible Navier-Stokes, most 
methods employ a fixed-point iteration solution strategy. Using fixed-point iteration provides a 
simple way to linearize the nonlinear convective term of the Navier-Stokes equations and also 
is relatively easy to implement. An alternative approach employs Newton-Raphson iterations to 
obtain the solution. As opposed to fixed-point iteration, for Newton-Raphson, one must construct 
the tangent matrix by taking derivatives of the weak form, including the nonlinear convective term. 
Although there are a number of drawbacks to using Newton-Raphson, the primary advantages 
include quadratic convergence and a framework on which to implement advanced algorithms such 
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as continuation methods [I]. The development of nomotopic or arc-length type methods shows 
potential for capturing high Reynolds flows otherwise unattainable with existing methods [2]. 

Obtaining numerical solutions to the Navier-Stokes equations intrinsically involves a number of 
challenges. The most formidable include instabilities in the pressure that arise due to the underlying 
Stokes problem, potential instabilities in the velocity for high Reynolds flows, and dealing the with 
nonlinearity of the convective term. In this work, we use a mixed finite element formulation that 
treats the velocity and pressure as independent variables. Although mixed formulations require 
more degrees of freedom and result in discrete systems that are indefinite (they have positive and 
negative eigenvalues), they pose the benefits of practicality and robustness in extreme cases [3]. 

The numerical instabilities that arise stem from the saddle point nature of the underlying Stokes 
problem. This feature is mathematically explained by the Ladyzhenskaya-Babuska-Brezzi (LBB) 
stability condition [J] . It is well known that the standard Galerkin formulation does not yield stable 
results, but a number of stabilized formulations have been developed to combat this deficiency. For 
a comprehensive overview of these methods see [5j El El [2] and the references therein. We generate 
both the fixed-point iteration and Newton-Raphson strategies in the context of a variational mul- 
tiscale framework. The variational multiscale concept [8] acknowledges a duality of scales inherent 
in the physics of the problem: a coarse or resolvable scale, and a fine-scale behavior that cannot 
be accurately captured by the mesh. To remedy instability, we model the fine-scale behavior by 
splitting the problem into two subproblems a coarse-scale subproblem and a fine-scale subproblem. 
We then construct our formulations in a way that models the fine-scale behavior using bubble 
functions. For the fixed-point iteration strategy, we solve the fine-scale subproblem in terms of the 
coarse-scale variables and substitute the result into the coarse-scale subproblem. For the consistent 
Newton-Raphson technique, we treat the two subproblems as separate nonlinear residuals and solve 
the system using an iterative update approach. 

The primary contributions of this work include the development of a consistent, stabilized, 
Newton-Raphson type formulation for incompressible Navier-Stokes and a direct comparison with a 
fixed-point iteration strategy. This work also lays the groundwork for advanced solution algorithms 
such as arc-length or homotopy methods that require origins in a consistent linearization. The 
work presented herein is the extension of work presented in [91 [TO] . 

In the following sections, we introduce the variational multiscale concept, followed by a formu- 
lation of the fixed-point iteration strategy. We then propose a Newton-Raphson solution technique 
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and derive the associated tangent stiffness matrix. The results presented focus on comparing the 
Newton-Raphson technique to fixed-point iteration for a number of example problems. For all of 
the example problems, we employ linear triangular elements. In [TT], we show that the consis- 
tent formulation is unstable for linear quadrilateral and hexahedral elements. We conclude with a 
summary of the main findings of this work. 



2. GOVERNING EQUATIONS 

Let ft be a bounded open domain, and T be its boundary, which is assumed to be piecewise 
smooth. Mathematically, T is defined as T := ft — ft, where ft is the closure of ft. Let the velocity 
vector field be denoted by v : ft — > M. nd , where "nd" is the number of spatial dimensions. Let 
the (kinematic) pressure field be denoted by p : ft — > K. As usual, T is divided into two parts, 
denoted by T v and r , such that n T* = and T v U T* = T. T v is the part of the boundary on 
which velocity is prescribed, and r is part of the boundary on which traction is prescribed. The 
governing equations for incompressible Navier-Stokes flow can be written as 

(1) v -Vv - 2vV 2 v + Vp = b in ft 

(2) V • v = in n 

(3) v = v p on 

(4) -pn + v(n ■ V)v = t n on r* 



where v is the velocity, p is the kinematic pressure (pressure divided by density), V is the gradient 
operator, V 2 is the Laplacian operator, b is the body force, v > is the kinematic viscosity, v p is 
the prescribed velocity vector field, t n is the prescribed traction, and n is the unit outward normal 
vector to V. Equation ([1]) represents the balance of linear momentum, and equation ([2]) represents 
the continuity equation for an incompressible continuum. Equations ([3]) and dH) are the Dirichlet 
and Neumann boundary conditions, respectively. 



2.1. Classical mixed formulation. First, let us define the relevant function spaces for the ve- 
locity, v(x) and the weighting function associated with velocity, which is denoted by w(x). They 
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are, respectively, defined as 

(5) V := {v | v G (H\n)) nd ,v = v p our"} 

(6) W:= {w|wG(ff 1 (fi)) ,!d ,w = Oonr} 

where is a standard Sobolev space [4]. In the classical mixed formulation, the function space 

for the pressure p(x) and its corresponding weighting function q(x) are given by 

(7) V:={p\peL 2 (n)} 

where L 2 (0) is the space of square-integrable functions on the domain Q. In the fixed-point iteration 
based variational multiscale method, the function space for p(x) and q(x) will be 

(8) V:={p\peH l {U)} 

For further details on function spaces refer to Brezzi and Fortin [3]. 

Remark 1. When Dirichlet boundary conditions are imposed everywhere on the boundary, that is 
r* = 0, the pressure can be determined only up to an arbitrary constant. In order to define the 
pressure field uniquely, it is common to prescribe the average value of pressure, 

(9) ( pdn = Po 

where po is arbitrarily chosen (and can be zero). Then, the appropriate function spaces for the 
pressure that should be used instead ofV( defined in equation (J7J) ,) is 

(10) V :={p\p£ L 2 (Q), [ p dU = 0} 

Jn 

Another way to define the pressure uniquely is to prescribe the value of the pressure at a point, 
which is computationally the most convenient. 

For convenience, we define the I? inner-product over a spatial domain, K, as 

(11) (a,b) K = I a b&K 

Jk 

The subscript K will be dropped if K is the whole of O, that is K = Vt. The classical mixed for- 
mulation (which is based on the Galerkin principle) for the incompressible Navier-Stokes equations 
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can be written as: Find v(x) £ V and p{x) € V such that 



(12) 



(to, v ■ Vw) + (Vio, 2i/Vo) - (V ■ io,p) = (to, b) + (to, h)r t V to € W 



(13) 



(g,V-«)=0 Vg£P 



To determine a numerical solution using the finite element method, one first chooses the ap- 
proximating finite element spaces, which (for a conforming formulation) will be finite dimensional 
subspaces of the underlying function spaces of the weak formulation. Let the finite element func- 
tion spaces for the velocity, the weighting function associated with the velocity, and the pressure 
be denoted by V h C V, W h C W, and V h C respectively. The finite element formulation of the 
classical mixed formulation then reads: Find v h (x) G V h and p (x) € V h such that 



For mixed formulations, the inclusions V C V, W C W, and P C V are themselves not sufficient 
to produce stable results, and additional conditions must be met by these finite element spaces to 
obtain meaningful numerical results. A systematic study of these types of conditions on function 
spaces to obtain stable numerical results is the main theme of mixed finite elements. One of the 
main conditions to be met is the LBB inf-sup stability condition. For further details, see [1J[2]. 

3. VARIATIONAL MULTISCALE FRAMEWORK 

To remedy the inadequacies of the standard Galerkin formulation as presented above, both 
the fixed-point iteration based and the Newton-Raphson based solution strategies incorporate a 
variational multiscale framework for stabilization. The variational multiscale concept, which stems 
from the pioneering work by Hughes [8], decomposes the underlying fields into coarse or resolvable 
scales and subgrid or unresolvable scales. This decomposition provides a systematic way to develop 
stabilized methods. 

3.1. Multiscale decomposition. Let us divide the domain f2 into N non-overlapping subdomains 
Q e (which in the finite element context will be elements) such that 



(14) (w h , v h ■ Vv h ) + (Vto h , 2uVv h ) - (V • w h ,p h ) = (w h , b) + (w h , h) Tt V w h G W 



(15) 



{q h ,V -v h ) =0V q h EP 



,h 



N 



(16) 




e=l 
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The boundary of the element Q e is denoted by T e . We decompose the velocity field v(x) into 
coarse-scale and fine-scale components, indicated as v(x) and v'(x), respectively. 

(17) v(x) = v(x) + v'(x) 

Likewise, we decompose the weighting function w(x) into coarse-scale w(x) and fine-scale w'(x) 
components. 

(18) w(x) = w(x) +w'(x) 

We further make an assumption that the fine-scale components vanish along each element boundary. 

(19) v\x) = w'{x) = on r e ; e = 1, . . . ,N 

Let V be the function space for the coarse-scale component of the velocity v, and VV be the function 
space for w; and are defined as 

(20) V:= V; W:=W 

where V and W are defined earlier in equation ([5]) and equation ([6]) respectively. Let V' be the 
function space for both the fine-scale component of the velocity v' and its corresponding weighting 
function w', and is defined as 

(21) V :={v\v£ (H 1 (Q e )) nd , v = Oonr e ,e = 1,...,N} 

The velocity field v(x) is now an element of the function space generated by the direct sum of 
V and V', denoted by V © V. Similarly the direct sum of VV and V', denoted by VV © V', is the 
function space for the field w(x). 

In theory, we could decompose the pressure field into coarse-scale and fine-scale components. 
However, for simplicity we assume that there are no fine-scale terms for the pressure p(x) and for 
its corresponding weighting function q(x). Hence, the function space for the fields p(x) and q(x) 
is V. 
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3.2. Two-level classical mixed formulation. The substitution of equations (|17p and (|18p into 
the classical mixed formulation, given by equations (|12p and (|13p . becomes the first point of depar- 
ture from the classical Galerkin formulation. 

(to + w', (v + v') ■ V(v + «')) + (V(tu + to'), 2i/V(« + «')) 

(22) _(v ■ (to + to'),p) = (to + to', 6) + (to + to', h) Tt 

(23) (q,V-(v + v'))=0 

Because the weighting functions to and to' are arbitrary, and because the functionals are linear in 
the weighting functions, we can write the above problem as two sub-problems. The coarse-scale 
problem can be written as: 

(to, (v + v') • V(t) + v')) + (Vto, 2vV{v + v 1 )) - (V • w,p) 

(24) = (w, b) + (w, h) Tt VtueW 

(25) (q, V • (v + v')) = Vq£V 
The fine-scale problem can be written as: 

(to', (v + v') ■ V(v + v')) + (Vto', 2z^V(^3 + «')) - (V • tu',p) 

(26) = (to', 6) + (to', fc) rt Vw'elV 

Remark 2. iVoie i/iai t/ie /ine scale problem is independent and uncoupled at the element level 
(defined over the sum of element interiors). Due to the assumption that the subgrid scale response 
vanishes on the element boundaries, (w,v + v') = (w,v) + (w,v'). 

3.3. Fine— scale interpolation and bubble functions. If one chooses a single bubble function 
for interpolating the fine-scale variables (similar to the MINI element [12]), then we have 

(27) v' = b e (3; w' = 6 e 7 

where b e is a bubble function, and (3 and -y are constant vectors. The gradients of the fine-scale 
velocity and weighting functions are 

(28) W = f3Vb eT ; Vto' = 7 V6 eT 

where V6 e is a dim x 1 vector of the derivatives of the bubble function. Standard bubble functions 
for several elements are provided in Table HJ 



Table 1 . Bubble Functions for Standard Finite Elements 



Element Bubble function 

66(1-6-6) 
666(1-6-6-6) 
(i-6 2 )(i-6 2 ) 
(i-6 2 )(i-e 2 2 )(i-6 2 ) 



4. FIXED-POINT ITERATION 

Equation ([1]) is nonlinear as engendered by the convection term v ■ Vv. One popular method for 
dealing with the nonlinearity of this term is to linearize the weak form using fixed-point iteration. 
For a derivation of the linearization using fixed-point iteration, see Appendix 18.21 The linearized 
forms of the fine and coarse-scale subproblems are written as: 
Linearized coarse-scale subproblem 

(w,v c ■ V(v + v')) + (w, (v + v') ■ Vv c ) + (Vw,2vV(v + v')) - (V • w,p) 

(29) =(w,b) + (w,h) Tt 

(30) (g,V- (v + v')) = 

Linearized fine-scale subproblem 

(w', v c ■ V(v + v')) + (w', (v + v') ■ Vv c ) + (V«)', 2vV(v + v')) - (V • w',p) 

(31) =(w',b) 

To derive our stabilization parameter, we eliminate the fine-scale variables by solving the fine-scale 
problem (equation (|3ip ) in terms of the coarse-scale variables. We then substitute the fine-scale 
solution into the coarse-scale problem (equation (|29p ) and solve the coarse-scale problem to obtain 
v(x) and p(x). Typically, the stabilization parameter is derived in a consistent manner by incorpo- 
rating the coarse-scale residual evaluated over the element. Examples of such formulations include 
the work of Masud and Khurram [13] for the Navier-Stokes equations and that of Nakshatrala et 
al [H| for nearly incompressible linear elasticity. The derivation, which is presented in Appendix 
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T3 

TET4 

Q4 

B8 



18.31 leads to a stabilized weak formulation of the following form: 



(w,v c • Vv) + (w,v- Vv c ) + (Vw,2uVv) - (V • w,p) + (q, V • «)+ 



( 



t> c • Vid + 2^V 2 w + Vg - (Vu c ) 1 w, r 



(Vjp + f c • Vu + v • Vt> c - 2v\7 2 v 



*>)) 



(32) 



(id, 6) + (w,h)r t 



Remark 3. Notice that we have introduced the laplacian of the weighting function and the velocity 
as a result of integrating by parts. Computing these terms requires implementing the third order 
tensor which represents the second derivative of the shape functions. 

Remark 4. Recall that these equations have been linearized about a certain velocity v* with change 
in velocity 5v. For simplicity, we dropped the superscript on v* and changed 5v to v c , so to recover 
the stabilized nonlinear form of the weak formulation we must "reverse" the linearization process. 
Given that even the stabilization parameter r contains v c and the complex expressions created 
from using integration by parts to isolate v' , recovering the nonlinear form of Equation\3S\ remains 
difficult if not impossible. 

Remark 5. Notice that the consistently linearized mixed variational multiscale form of the Navier 
Stokes equations is not self- adjoint as shown by the negative sign of the fourth term in the weighting 
slot of the stabilization term. 

The fixed-point iteration based weak form (equation (|32|) ) is useful for low and moderate Reynolds 
flows, but the rate of convergence may be lower and the take a greater number of iterations [2]. In 
order to utilize homotopy methods (similar to arc-length), which greatly increase convergence for 
high Reynolds flows, one must posses the nonlinear weak formulation. 



Instead of linearizing the convective term and solving the coarse and fine scale subproblems by 
substitution, as in the fixed point iteration strategy, for the consistent Newton-Raphson method, 
we treat equations ([24")) - ([26]) as global residuals and obtain the solution using an iterative update 
equation. 



5. CONSISTENT NEWTON-RAPHSON SOLUTION STRATEGY 
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5.1. Scalar residual. The scalar residual equations may be written as 



r c (v;v',p) :=(w, (v + v') ■ V(v + v')) + (Vw,2vV(v + v')) 

(33) - (V -w,p) - (w,b) - (w,h) Ti 

(34) r p (v;v') := - (q, V • (v + v')) 

(35) r f (v;v',p) :={w', (v + v') ■ V(« + v')) + (Vio', 2i/ / V(« + «')) ~ ( v ■ "Ap) ~ C™') b ) 
where the subscripts 'c', 'p', and 'f stand for coarse, pressure, and fine. 

5.2. Vector residual. We again choose a single bubble function for interpolating the fine-scale 
variables. The solution, v, and its weighting function, w, may be expressed in terms of the nodal 
values v and w as 

(36) v = % T N T ; w = -0) T N T 

where AT is a row vector of shape functions for each node. Substituting equations (|27|) and (|36p 
into equations (|33p - (|35p and noting the arbitrariness of w and 7, we construct vector residuals, R, 
that are the sum contributions of the vector residuals at the element level, R e , given as 

R e c (v; v',p) := [ (N T I)((v + v') -V)(v + v') dft 
Jn 

+ 2v I ((DN)J~ l J)vec[Vv + Vi/] dQ 
Jn 

(37) - f vec[J~ T (DN) T ]pdn- [ (N T I)b dQ, — [ (N T I)h dT 

Jn Jn Jv t 

(38) Ri(v;v'):=- [ N T V ■ (v + v') 

Jn 

R e Jv;v',p) := [ b e I((v + v 1 ) ■ V)(v + v') 
Jn e 

+ 2v I (V6 eT © J)vec[Vu + W] dft 
Jn e 

(39) -I (Vb eT Q I)vec[I]p dtt - f b e Ibdft 



dft 

dn 
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where DN represents a matrix of the first derivatives of the element shape functions, which is 
defined for a triangular element as 



(40) 



DN 



dNi dNi 
9£i d& 



8N 3 dN 3 
L 9£i d& 



J is the element jacobian matrix, vec[-] is an operation that represents a matrix with a vector, 
is the Kronecker product [15] (see Appendix 18. 

5.3. Extension to transient flows. To extend the formulation presented above to transient flows, 
we need only to incorporate the time derivative of the velocity, v, in the balance of momentum, 
equation ([I]). This leads to the addition of the following terms, (w,v) and (w',v) to the residual 
equations (|33p and (|35p . respectively. This in turn adds the following terms, J n (N T I)v dQ and 
f n b e Iv dQ, to the vector residuals, equations ([37|) and ([39]) . respectively. 

To integrate in time, from step n to step n + 1, we use the backward Euler method by which we 
express the acceleration, v, as the average of the velocity over time. 

v n+ i - v n 



(41) 



V n +1 



At 



Substituting equation ([41]) into the additional acceleration terms in ([37|) and ([39]) . we have 
(42) I (N T I)i, n+1 dQ= f ^- t (N T I)v n+1 dn - f -^(N T I)v n dO 



(43) 



6 e /^ dn = / ^iwn+i dn - / -^ t iv n dn 



where v n is the converged velocity from the previous time step. 

5.4. Tangent matrix. Using a Newton-Raphson type approach, we obtain the solution in an iter- 
ative fashion using the following update equation until the residual is under a prescribed tolerance, 



(44) 



P 



i+i 



jH+l = V H + Av f 



where the updates at each iteration, % are calculated from the following system of equations 



(45) 





DR C 


DRc 


Dv 


Dp 


Dv' 


DRp 


DRp 


DRp 


Dv 


Dp 


Dv' 


DR f 


DR f 


DR f 


Dv 


Dp 


Dv' 



Av 






Ap 


H 


Rp 


Av' t 




k R f 



The element matrices, j^j, which are assembled to form the consistent tangent matrix, are derived 
in Appendix 18.41 

6. NUMERICAL RESULTS 

In this section, we present the results from several numerical experiments that directly compare 
the Newton-Raphson solution approach to fixed-point iteration. 

6.1. Body force driven cavity. We first perform a mesh refinement study to observe the con- 
vergence rates for the Newton-Raphson and fixed-point iteration strategies for a smooth problem 
using linear elements as the element length decreases. The problem consists of a body force driven 
cavity for which the solution may easily be obtained analytically. The geometry and a typical mesh, 
shown in Figure Q3 consists of a unit domain with boundaries that prohibit flow in the tangential 
and normal directions (v x = v y = 0.0). The prescribed constant body force, b, in components, is 
given as 

b x = (12 - 24y)x 4 + (-24 + 48y)x 3 + (-48y + 72y 2 - 48y 3 + 12)x 2 
+ (-2 + 24y - 72y 2 + 48y 3 )x + 1 - 4y + I2y 2 - 8y 3 

+ (4xy - 12xy 2 + 8xy 3 - 12x 2 y + 36x 2 y 2 - 24x 2 y 3 + 8x 3 y - 2Ax 3 y 2 + 16x 3 y 3 )* 
(x 2 (l-x) 2 (2y-Qy 2 + 4y 3 )) 

+ (2x 2 - 12x 2 y + 12x 2 y 2 - Ax 3 + 24x 3 y - 24x 3 y 2 + 2x 4 - 12x 4 y + 12x 4 y 2 )* 

(46) (-y 2 (l-y) 2 (2x-6x 2 + Ax 3 )) 

by = (8- A8y + 48y 2 )x 3 + (-12 + 72y - 72y 2 )x 2 

+ (4 - 24y + 48y 2 - 48y 3 + 24y 4 )x - I2y 2 + 24y 3 - 12y A 

+ {-2y 2 + I2y 2 x - 12y 2 x 2 + 4y 3 - 24y 3 x + 24y 3 x 2 - 2y 4 + 12y 4 x - 12y 4 x 2 )* 
(x 2 (l-x) 2 (2y-6y 2 +4y 3 )) 

+ (-4yx + 12yx 2 - 8yx 3 + 12y 2 x - 36y 2 x 2 + 24y 2 x 3 - 8y 3 x + 24y 3 x 2 - 16y 3 x 3 )* 

(47) (-y 2 (l-y) 2 (2x-6x 2 + 4x 3 )) 
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For a unit viscosity, the exact solution is 

v x = x 2 (l-x) 2 (2y-6y 2 + 4y 3 ) 
Vy = - y 2 (l - y) 2 (2x - 6x 2 + Ax 3 ) 
(48) p = x(l - x) 

The results, for triangular elements of edge length h with a prescribed body force and known 
solution, are shown in Figures [2] and [3l In Figure [2] the values of the pressure and velocity are 
shown as measured along the x-axis at y = 0.5. From Figure [3] we see that the i? 1 -semi norm of 
the pressure is 1.0 and the L 2 norm of the velocity is 2.0, respectively for both methods. We also 
see that the error in the velocity is slightly lower for the Newton-Raphson strategy. 

In Figured] we show a comparison between the residual convergence for the consistent Newton- 
Schur approach and the fixed-point iteration method for the smooth body force problem at Re = 
400. The solution converges in three steps for the Newton-Raphson approach, whereas after ten 
iterations, the fixed-point method has not yet reached the tolerance. For the higher Reynolds case, 
Re = 5, 000, the Newton-Raphson strategy again takes three iterations to converge, but the fixed- 
point iteration method diverges. In this case, one would have to employ a continuation technique 
to obtain the solution using the fixed-point iteration method. 

6.2. Lid-driven cavity. The lid-driven cavity has served as a benchmark problem for numerical 
methods and has been analyzed by a number of authors (for examples, see |16t [T7\ \TE[ I19|). The 
problem description, boundary conditions, and a typical mesh are shown in Figure For a unit 
length and height domain, a unit velocity in the horizontal direction is prescribed across the top 
boundary, while the sides and bottom are treated as solid surfaces with no flow in the horizontal 
or vertical direction. To resolve the boundary conditions at the corners they are treated as part of 
the vertical walls. Therefor, the prescribed velocity approaches zero at the corners. Plots of the 
stream traces over the velocity contours are shown for Re = 400 and Re = 5, 000 in Figures [7] and [8] 
respectively for the Newton-Raphson approach. Similar plots for the fixed-point iteration strategy 
have been presented in [13]. Results for both methods are compared to other published results in 
Figure [6l Figure [6] shows the velocity magnitude in the x-direction as measured along the y-axis at 
x = 0.5. The results closely match those presented in [16]. For the Re = 5, 000 case, aside from the 
central vortex, the lower vortices and the upper left vortex are captured for a mesh of triangular 

13 



elements with edge length 0.025 using the Newton-Raphson based approach. We were unable to 
obtain a converged solution to the lid-driven cavity problem using the fixed-point iteration scheme 
at Re = 5, 000. 

For the lid-driven cavity, we observe quadratic convergence as expected for the consistent Newton- 
Raphson approach. The residual is listed for each iteration in Table [2j Above Re = 5, 000, we 
were unable to obtain a converged solution for the Newton-Raphson based approach. Even when 
a continuation method is employed by solving the problem at low Re and slowly increasing Re 
by a factor of 1.1, using the previous converged solution as the starting point, the method does 
not converge. There seems to exist a well defined Re above which the consistent method does not 
converge. The lack of convergence for consistent methods at very high Re is pointed out in [2]. 

Table 2. Quadratic Convergence for Lid-Driven Cavity Problem at Re = 400 



Iteration 


Residual 


1 


1.429492E+04 


2 


8.608793E+03 


3 


4.716364E+03 


4 


2.104800E+03 


5 


6.737382E+02 


6 


5.754620E+01 


7 


4.170450E-01 


8 


5.516418E-04 


9 


4.784608E-08 


10 


1.784712E-11 



6.3. Backward facing step. Another well-known benchmark problem that posses a corner sin- 
gularity is the backward facing step. The domain, boundary conditions, and mesh are shown in 
Figure [H A parabolic velocity is prescribed over the inflow face with a maximum value of 1.0 at 
the center of the opening. The results for the consistent Newton-Raphson method are presented 
in Figure [TOl which show streamtraces over the velocity contours and the pressure contours respec- 
tively for Re = 150. Results for the fixed-point iteration method are presented in [13]. In order to 
obtain the solution for both the fixed-point iteration and consistent Newton-Raphson approaches, 
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a continuation process must be employed. To compare the convergence of the residual for both 
cases, we present the results for the first converged Reynolds number in the continuation. Figure [TT1 
shows the convergence of the residual at Re = 15. Again, the Newton-Raphson method converges 
in much fewer iterations. 

6.4. Consistent Newton-Raphson and the transition from steady to periodic flow. 

Whereas the convergence for the Newton-Raphson based approach is much faster than the fixed- 
point iteration based method, the Newton-Raphson based approach exhibits vulnerability in the 
region of turning or bifurcation points. In particular, the consistent method does not converge when 
the flow past a bluff object transitions from steady to periodic vortex shedding. This transition 
represents an instance of Hopf bifurcation [20] . 

The domain and boundary conditions for two problems that engender this behavior are shown 
in Figures [12] and [T3] for flow past a circular and square cylinder respectively. The computational 
meshes used for this problem are also shown in Figures [12] and [ill Figures [13] and [15] show the last 
converged iterations for the consistent Newton-Raphson based approach immediately before the 
method fails to converge. Notice that for both cases, periodic vortex shedding has not yet started. 
The residual for each case is shown in Tables [3] and [4] 

The consistent Newton-Raphson method does not fail for all transient vortex shedding related 
problems, for example, the jet flow through an orifice problem described in Figure [16] The com- 
putational mesh is also shown in Figure [16] Figure [17] shows the results using the consistent 
Newton-Raphson approach for a jet flowing through a 1/16 unit opening in an infinite wall for vis- 
cosity, v = 0.005, and maximum flow speed at the orifice of v = (1.0,0). The results closely match 
those reported in [5] for the same Reynolds number. Again quadratic convergence is obtained for 
the Newton-Raphson technique for an initial guess of v = 0. For the jet flow problem, there is no 
transition from steady to periodic flow. 

7. CONCLUSIONS 

We have presented a comparison between a fixed-point iteration and consistent Newton-Raphson 
solution strategy for a variational multiscale formulations for incompressible Navier-Stokes. The 
Newton-Raphson approach shows faster convergence for several numerical test problems and in 
many cases converges quadratically in the vicinity of the solution. We have shown that in some 
instances, the Newton-Raphson approach converges while the fixed-point iteration strategy does 
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Table 3. Convergence for Flow Past a Circular Cylinder at Re = 75 



Circular Cylinder 


Iteration 


Prior to Hopf Bifurcation 
Residual 


At Hopf Bifurcation 
Residual 


1 

2 
3 
4 
5 
6 
7 


7.338022E+01 
5.783447E-02 
4.012763E-02 
4.246893E-07 


6.066543E+01 

2.542188E-01 

2.447873E-01 

2.451840E-01 

2.451840E-01 

2.451840E-01 


Table 4. Conver; 


j;ence for Flow Past a Squar 


e Cylinder at Re = 75 


Square Cylinder 


Iteration 


Prior to Hopf Bifurcation 
Residual 


At Hopf Bifurcation 
Residual 


1 

2 
3 
4 
5 
6 
7 


2.526981E+01 
1.561594E-02 
5.892473E-03 
4.721279E-08 


1.460652E+00 

7.228785E-02 

7.196581E-02 

7.195860E-02 

7.195860E-02 

7.195860E-02 



not. We have also shown that for some problems, the Newton-Raphson strategy is unable to pass 
through the transition from steady to periodic flow without additional sophistication being built 
into the algorithm. In summary, this work points to the need for advanced solution algorithms able 
to maintain the convergence qualities of the Newton-Raphson technique with increased robustness 
for high Reynolds flows and flows with distinct bifurcation points. 



8. APPENDIX 
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8.1. Notation and definitions. Consider an n x m matrix A and ap x q matrix B 



CL n I ... flj; 



'1,1 



3p,l 



The Kronecker product of these matrices is an np x mg matrix, and is defined as 



AQB :-- 



o>i,lB . . . a^ m B 



a n iB . . . a nm B 



The vec[-] operator is defined as 



vec[A] :-- 



8.2. Fixed point linearization of the convective term. We linearize the convective term, 
v ■ Vu in the weak form about a velocity v* with a change in velocity 6v as follows. 



r(v) = (w, v ■ Vv) 



d_ 

dT 



r(v* + edv) 



d r 
=o de 



€=0 



(49) 



(w, (v* + e5v) ■ V(v* + e8v)) 
(w, 5v ■ V(v* + edv)) + (w, (v* + e5v) ■ V5v) 
(w, 5v ■ Vv*) + (w, v* ■ V5v) 



e=0 



To simplify the formulation, we drop the superscript on v* and refer to the change in velocity 5v 
as v c . The linearized form of the second term on the left side of equation f|12 j) is now expressed as 



(50) 



(w, v ■ Vv) = (w,v c ■ Vv) + (w,v ■ Vv c ) 
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8.3. Derivation of the stabilization term for the fixed-point iteration strategy. Due to 

the linearity of the solution and using integration by parts on the pressure term, we re-arrange 
equation (f3"Tj) as 

(w', v c ■ Vv') + (w', v' ■ Vv c ) + (Vw', 2vVv') = (w', b f ) - (w', Vp) 
(51) ~(w', v) - (w', v c ■ Vv) - (w', v ■ Vv c ) - (Vv/, 2uVv) 

Performing integration by parts on the last term in the previous equation and collecting terms we 
reduce this expression to 

(52) 

(«/, v c ■ Vv') + (w', v' ■ Vv c ) + (Vw', 2uVv') = (w', b f - Vp - v - v c ■ Vv - v ■ Vv c + 2uV 2 v) 
Making the substitution r = bf — Vp — v — v c ■ Vv — v ■ Vv c + 2vV 2 v, we have 

(53) (w', v c ■ Vv') + (w', v' ■ Vv c ) + (Vw', 2uVv') = (w', f) 

We now solve equation (|53p analytically using bubble functions. We assume that v' and w' are 
represented by bubble functions over Q'^, as in equation (|27p . and move the constant coefficients 
outside of the integrals to obtain 

(54) 7 T ^(b e v c - Vb e + v\Vb e \ 2 )I+(b e ) 2 Vv c + vVb e ®Vb e d{?J (3 = 1 T (b e ,f) 

where Vv c is the gradient of the known velocity about which the equations have been linearized 
and Vb e is a dim x 1 vector of the derivatives of the bubble functions. Because the coefficients 7 
are arbitrary, we have a system of equations that we can solve for f3 in terms of f , 

(55) f3 = A^R 
where A and R are defined as follows 

(56) A= f (b e v c ■ Vb e + v\Vb e \ 2 )l + (b e ) 2 Vv c + vVb e ®Vb e &SI 

(57) R= / b e rdn 

Jne f 

we now have an expression for the fine-scale velocity which we can substitute into the coarse-scale 
subproblem. 

(58) v'(x) = b e (3 = FA-iR 
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In the finite element setting, as the mesh becomes adequately refined, r will essentially be constant 
over an element. We can then move r to the outside of the integral in equation (|57p and rewrite 
our approximation for the fine-scale velocity as 



v'(x) = b e b e dfiA^f = rr 



(59) 



where the stabilization parameter r (in this case a second-order tensor) is defined as 

T = b e [ b e dnA- 1 



f 



(60) 



3 [ b e dn [ (b e v c -Vb e + v\Vb e \ 2 )I+(b e ) 2 Vv c + vVb e ®Vb e dQ 
Jn% Jn* 



Using integration by parts, the linearity of the solution field, and the identity, (a, Ab) = (A T a, b), 
we substitute equation (|59p into the coarse scale subproblem to get 

(w,v) + (w,v c • V«) + (w,v- Vv c ) + {Vw,2vVv) - (V • w,p) + (q, V • v) + 
(v c ■ Vii) + 2vV 2 w + Vq- (Vv c ) T w, r(Vp + v + v c ■ Vv + v ■ Vv c - 2vV 2 v - b)) 

(61) ={w,b) + (w,h) rt 

8.4. Derivation of the terms for the consistent tangent matrix. 

5^ • Sv = ! (N T I)Sv dn + I (N T I)((Sv ■ V)(« + v') + (v + «') • V5u) dft 

(62) + 2u f [ ({DN)J- l Ql)vec[VSv]dn 

Jn 



(63) 



5^ • Sv = - I N T V ■ Sv dfi 
Di> Jn 



(64) 
(65) 
(66) 



j ■ Sv' = / (N T I)((Sv' ■ V)(w + «') + (« + t/) • Vtfw') dft 
Jn 

+ 2v f [ ((DN)J~ 1 I)vec[V5v'] dU 
Jn 



BR 



p • Sv' 



Dv' 



N T V ■ Sv' dtt 



BR C 
Bp 



Sp = - [ vec[J- T (DN) T ]Spdn 
Jn 
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DR 



L . 5v = I b e I5T> dtt + I b e I{8v • V)(« + v') + ((« + «') • V)«5« dtt 



D.R, 



(67) + 2v f [ (Vb eT Ql)vec[V8v]dtt 

Jn e 

if r 

-f ■ Sv' = / & e J(<fo' • V)(« + «') + ((« + »') • V)«Jw' dO 
; Jn e 

(68) + 2i// / (Vb eT Ql)vec[V8v']dtt 

(69) • 8p = - I (Vb eT I)vec[I]Sp 

Dp Jn e 

In the finite element setting, we introduce a finite element discretization of the solution for the 
pressure and velocity as follows 

(70) v = t) T N T ; v' = b e {£)(3 ; p = Np 
Making use of the following equivalences: 

(71) m.„ m %.^..m.„ m w. Vi m.* m w. 9 

we derive the the discretized form of the tangent stiffness matrix as follows, 

2^ • vec[^ T l = ( (N T N I)^-vec[8f> T ] dtt + ( (N T I)V(v + v')(N J)vec[<^ T ] dtt 
Du J n At J n 

+ f {N T (v + v') T J~ T {DN) T I)vec[8t> T ] dtt 
Jn 

(72) + 2v f [ (( J DAT)J- 1 J- T ( J DAT) T 0/)vec[^ T ] dtt 

Jn 

(73) • vec[^ T ] = - / N T vec(J" T (DN) T ) T vec[5v T ] dQ 
Dv Jn 

2^ • 8(3 = [ (N T I)b e V(v + «') 5/3 dO + f (N T I)Vb eT (v + v')I 8(3 dtt 

(74) + 2i/y / {(DN)J- 1 Vb e J) 5/3 dO 

Jn 



(75) ^§ • 5/3 = - / N T Vb eT 8(3 dtt 

Dp Jn 

(76) P^p.<^ = - / vec[J _T (i?iV) T ]iV(Jpdn 

Dp Jn 
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. vec[5t> T ] = I b e I(N 7)^-vec[5^ T ] dft + I b e IV(v + v')(N 7)vec[5f> T ] dft 
+ f b e I((v + v') T J- T (DN) T 7)vec[5l> T ] dft 

(77) + 2v f I (Vb eT J~ T (DN) T Ql)vec[V5v] dft 

(78) — / • 5/3 = / b e I(b e V(v + w') + (V6 e • (« + v'))I)5f3 dft + 2i// / (V6 eT W 7)5/3 dft 

(79) —L.8p = - / (V6 eT 7)vec[7]iV5pdft 

Dp 7q e 

D.R e 

Due to the nontrivial nature of the derivation of -j^- ■ 5(3, the derivation is presented as follows: 

^■5(3= [ b e I(b e 5(3 ■ V)(« + «') + ((« + «') • V)(6 e 5/3) dft 
Dp Vn e 

+ 2i// / (W T 7) (V6 e 7)5/3 dft 
= f b e I(b e 5(3 ■ V)(« + «') + (5/3 V6 e ) • (« + i/) dft 
+ 2v f [ (Vb eT 7)(V6 e 7)5/3 dft 

= [ b e I(b e V(v + v') + (V6 e • (v + v'))I)5/3 dtt 

Jn e 

(80) + 2i// / (Vb eT Vb e 7)5/3 dft 

Noting the arbitrariness of vec[5?) T ], 5/3, and dp, we have expressions for the tangent stiffness matrix 
components in the finite element setting: 

^ = J (N T N QI)±- + (N T Q I)V(v + v')(N 7) + (N T (v + v') T J~ T (DN) T 7) 

(81) + 2u f {{DN)J- l J- T {DN) T 7) dft 

(82) = / N T vec[J- T (DN) T ] T dft 

D^ 7n 

(83) = / (iV T I)b e V(v + «') + {N T Vb eT {v + v') 7)+ 2v f {{DN)J~ l Vb e 7) dft 
Di> 



(84) -^-f = / JV V6 dft 



n 
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(85) ^ = f -vec[J- T {DN) T ]N dtt 

Dp Jn 



[ b e I(N !)-*- + b e IV(v + v')(N I) + 6 e I((-u + v') T J~ T (DN) T /) 



Dv 



(86) + 2u f (Vb el J 1 (DN) /) 

(87) 1>J~ = J fee2v (^ + ^) + (V6 eT (^ + ^))^^ + 2^(V6 eT V6 e 0/)dO 



: / -(X7b eT Q I)vec[I]N dn 



Dp 
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Figure 1. Body force driven cavity: problem statement and boundary conditions 
(left) computational mesh (right). 




FIGURE 2. Body force driven cavity: pressure and velocity in the y-direction as 
measured along the x-axis at y = 0.5. 
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FIGURE 3. Body force driven cavity: H 1 -semi norm of the pressure vs. element size 
for Re = 400 (left) 1? norm of the velocity vs. element size (right). 
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FIGURE 4. Body force driven cavity: convergence of the residual for Re = 400 (left) 
Re = 5,000 (right). 
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FIGURE 5. Lid driven cavity: problem statement and boundary conditions (left) 
computational mesh (right). The non- leaky cavity approach is used here which 
resolves the discontinuity at the upper two corners of the domain by assuming that 
the corners belong to the vertical walls. 
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Figure 6. Lid driven cavity: velocity in the x-direction as measured along the 
y-axis at x = 0.5. 
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Figure 8. Lid driven cavity: velocity streamtraces over velocity magnitude con- 
tours for Re = 5, 000 using the consistent Newton-Raphson approach. 



27 




No-slip wall 

////// 



0,0 



1,0 f- 



? 7 7. / J 
No-slip wall 

3.0- 



1 



il 



Figure 9. Backward facing step: problem statement and boundary conditions (top) 
computational mesh (bottom). 




Figure 10. Backward facing step: using the consistent Newton-Raphson method 
for Re = 150, velocity streamtraces over velocity magnitude contours (top) pressure 
contours (bottom). 
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Figure 11. Backward facing step: convergence of the residual for Re = 15. 
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Figure 12. Flow past a circular cylinder: problem statement and boundary con- 
ditions (top) computational mesh (bottom). 



Correspondence to: Daniel Z. Turner, Department of Civil and Environmental Engineering, 2103 
Newmark Laboratory, University of Illinois at Urbana-Champaign, Urbana, Illinois - 61801. 
E-mail address: dzturnel@illinois.edu 

29 



Figure 13. Flow past a circular cylinder: streamtraces shown over pressure con- 
tours for Re = 75. 
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Figure 14. Flow past square cylinder: problem statement and boundary conditions 
(top) computational mesh (bottom). 
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Figure 15. Flow past a square cylinder: streamtraces shown over pressure contours 
for Re = 75. 
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Figure 16. Jet flow: problem statement and boundary conditions (left) computa- 
tional mesh (right). 
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